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ABSTRACT 

In a previous work, we have shown that the formation of the Fermi bubbles can be due to the 
interaction between winds launched from the hot accretion flow in Sgr A* and the interstellar medium 
(ISM). In that work, we focus only on the morphology. In this paper we continue our study by 
calculating the gamma-ray radiation. Some cosmic ray protons (CRp) and electrons must be contained 
in the winds, which arc likely formed by physical processes such as magnetic reconnection. We have 
performed MHD simulations to study the spatial distribution of CRp, considering the advection and 
diffusion of CRp in the presence of magnetic field. We find that a permeated zone is formed just outside 
of the contact discontinuity between winds and ISM, where the collisions between CRp and thermal 
nuclei mainly occur. The decay of neutral pions generated in the collisions, combined with the inverse 
Compton scattering of background soft photons by the secondary leptons generated in the collisions 
and primary CR electrons can well explain the observed gamma-ray spectral energy distribution. 
Other features such as the uniform surface brightness along the latitude and the boundary width of 
the bubbles are also explained. The advantage of this “accretion wind” model is that the adopted 
wind properties come from the detailed small scale MHD numerical simulation of accretion flows and 
the value of mass accretion rate has independent observational evidences. The success of the model 
suggests that we may seriously consider the possibility that cavities and bubbles observed in other 
contexts such as galaxy clusters may be formed by winds rather than jets. 

Subject headings: accretion, accretion disks-black hole physics-galaxies: active galaxies-galaxies: jets- 
Galaxies: nucleus 


1. INTRODUCTION 

Two giant gamma-ray bubbles above and below the 
Galactic plane were discovered by the Fermi Gamma- 
ray Space Telescope (Su et al. 2010; Dobler et al. 2010). 
Subsequently, additional observations have been per¬ 
formed and abundant data has been accumulated (Su 
& Finkbeiner 2012; Hooper & Slatyer 2013; Yang et al. 
2014; Ackermann et al. 2014). The main observational 
results are summarized as follows. 

• The bubbles extend to ~ 50° above and below the 
Galactic plane, and the width is ~ 40° in longitude. 

• The surface brightness is roughly uniform, with 
a significantly enhanced substructure in the south 
bubble, which is called the “cocoon”. 

• The boundary of the bubbles is sharp, with a width 
of about 3°. 

• The gamma-ray spectrum is uniform and hard. 
The spectral index in dN/dE oc U _T is 7 « 
1.9 ±0.2. The spectrum energy distribution (SED) 
shows a cut-off at ~ 100 GeV in high energy band, 
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and a significant cut-off below 1 GeV (Su et al. 
2010), although the latter is found to be less obvi¬ 
ous in a more recent study (Ackermann et al. 2014). 

• The total gamma-ray luminosity of both bubbles is 
4.4^4 x 10 37 erg s' 1 2 3 4 in 0.1 500 GeV band. 

Several models have been proposed to explain the for¬ 
mation of the Fermi bubbles. They usually invoke the 
interaction between the winds or jet and the interstellar 
medium to explain the morphology of the bubbles. These 
include: 1) winds radiatively driven from a quasar phase 
accretion disk in Sgr A* (Zubovas et al. 2011; Zubovas 
& Nayakshin 2012); 2) winds magnetically driven from a 
hot accretion flow in Sgr A* (Mou et al. 2014, hereafter 
paper I); 3) a jet launched from an accretion flow (Guo & 
Mathews 2012; Guo et al. 2012; Yang et al. 2012, 2013); 
4) star formation winds in the Galactic center (Crocker 
& Aharonian 2011; Crocker 2012; Crocker et al. 2014a,b; 
Carretti et al. 2013); 5) winds produced by the periodic 
star capture processes in Sgr A* (Cheng et al. 2011, 2014, 
2015). 

In terms of the origin of gamma-ray photons of the 
Fermi bubbles, theoretical models can be divided into 
two categories: the “leptonic” and “hadronic” ones. In 
the “leptonic” scenario, the gamma-ray photons come 
from the inverse Compton scattering (IC) of soft pho¬ 
tons (interstellar radiation field and cosmic microwave 
background) by relativistic electrons, which often called 
cosmic ray electrons (CRe). The origin of CRe is differ¬ 
ent in different models. They can come from: 1 )Fermi 
first-order acceleration in shock front formed in the pe¬ 
riodic star capture processes by Sgr A* (Cheng et al. 
2011, 2014, 2015); 2) the Fermi second-order accelera- 
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tion through stochastic scattering by plasma instabilities 
(Mertsch & Sarkar 2011); 3) directly carried by the jet 
(Guo & Mathews 2012; Guo et al. 2012; Yang et al. 2012, 
2013; Barkov & Bosch-Ramon 2014) or winds driven by 
the past star formation (Carretti et al. 2013). For the 
leptonic scenario, the cooling timescale of CRe with en¬ 
ergy of a few hundred GeV is only around 1 million years 
(Su et al. 2010). Hence, if the CRe come from galactic 
center (GC), i.e., carried by jet or wind from the accre¬ 
tion flow of Sgr A*, the age of the Fermi bubbles would 
be constrained to be less than 1-2 million years (Guo & 
Mathews 2012; Guo et al. 2012; Yang et al. 2012, 2013). 
Consequently, the power of the jet or wind must be rela¬ 
tively high so as to easily push the ISM away. This results 
in a high temperature of a few keV in the shocked ISM. 
Such a temperature is however higher than that detected 
in recent X-ray observations (Kataoka et al. 2013; Tahara 
et al. 2015). 

In the “hadronic” model, the gamma-ray photons are 
produced by the decay of neutral pions produced by the 
inelastic collisions between cosmic ray protons (CRp) 
and thermal nuclei in the ISM (pp collisions). The 
cooling timescale of CRp through pp collisions is t pp > 
3.5 x 10 9 yr (0.01 cm -3 ^ 1 ). Hence this model does not 
require a very high kinetic power from the GC. Another 
advantage is that, this model can naturally explain the 
drop of gamma-ray spectrum close to 1 GeV (Crocker & 
Aharonian 2011). 

There are many possible origins for the CRp. They 
may be: 1 ) injected in the wind (outflows) driven by the 
star formation in the GC, and have accumulated for an 
extremely long time (e.g., a few 10 8 — 10 9 years; Crocker 
& Aharonian 2011; Crocker 2012; Crocker et al. 2014a,b). 
2 ) accelerated in the accretion flow by processes like weak 
shocks and magnetic reconnection, and then carried into 
the bubbles by winds of accretion flow; 3) accelerated 
in and injected with the jet; 4) accelerated by the mul¬ 
tiple shocks or turbulence inside the Fermi bubbles, or 
the strong shock associated with the expansion of the 
bubbles (e.g., Fujita et al. 2013, 2014); 5) some diffu¬ 
sive injection of Galactic CRp and injection of CRp due 
to the expansion of the bubbles in the background ISM 
which is already full of CRp (Thoudam 2013). 

In paper I, we mainly focus on how to form the bub¬ 
bles, i.e., interpreting their morphology. By performing 
hydrodynamic simulations, we have shown that the bub¬ 
bles can be successfully produced by the interaction be¬ 
tween the winds launched from the hot accretion flow 
around Sgr A* and the ISM. The parameters of winds 
adopted in the model are taken from small scale MHD 
numerical simulations of accretion flows of Yuan et al. 
(2012). The required wind kinetic power, which is the 
main parameter of the model, is also in good consistency 
with that obtained from independent observational con¬ 
strains (see the review in Totani 2006). In addition, the 
model has also successfully explained the recent observa¬ 
tional results of ROSAT X-ray structure (Kataoka et al. 
2013). Most recently, the high-resolution X-ray observa¬ 
tion to the absorption line, combined with the ultraviolet 
observations, obtain the bulk motion velocity and tem¬ 
perature of the gas around the edge of the Fermi bubbles 
(Fang & Jiang 2014). Both of them are consistent with 
our simulation result. 


In the present work, we investigate the gamma-ray ra¬ 
diation of the Fermi bubbles. We adopt the “hadronic” 
scenario, i.e., the gamma-ray radiation comes from the 
pp collisions. For the origin of CRp, we consider the sec¬ 
ond mechanism mentioned above, i.e., CR protons are 
accelerated in the accretion flow (including the corona) 
and injected into the bubbles together with the winds. 
CR electrons must also be produced together with CRp. 
We will include their contribution when calculating the 
gamma-ray radiation. However, since the energy density 
of CRe is much lower than CRp, we can neglect CRe in 
numerical simulations (see §4.2). The main task is then 
to calculate the distribution of CRp. The motion of each 
CR particle is controlled by the magnetic field, so we 
should include magnetic field in our simulation. On the 
one hand, due to the scattering of magnetic irregularities 
in the winds, CRs are trapped in the plasma, i.e., they 
are “advected” by the wind gas. On the other hand, CRp 
can also diffuse with respect to the gas as they scatter off 
magnetic inhomogeneities. This will have an important 
effect on the CR distribution thus should be included in 
our simulation. In fact, as shown by Yang et al. (2012) 
and our present work, the anisotropic diffusion due to the 
configuration of the magnetic field lines is crucial to the 
interpretation of the sharpness of the observed surface 
brightness of the bubbles. To simulate the distribution 
of CRs together with the wind plasma, following Yang et 
al. ( 2012 ), we adopt a simplified approach, i.e., we treat 
CRs as a second fluid and solve directly for the evolution 
of CR pressure. 

The rest of the paper is organized as follows. In §2 we 
briefly review our accretion wind model. In §3, we in¬ 
troduce the details of our simulation setup. The results 
of the simulation and the calculation of gamma-ray ra¬ 
diation will be presented in §4. We then summarize and 
discuss our results in §5 and § 6 , respectively. 

2. MODELS OF ACCRETION WIND AND COSMIC RAYS 

2 . 1 . The Past Activity of Sgr A* and the Accretion 

Flow 

Sgr A* is extremely dim currently. However, many ob¬ 
servational evidences have been accumulated and shown 
that Sgr A* very likely was much more active during the 
past millions of years (see reviews by Totani 2006 and 
Bland-Hawthorn et al. 2013). The mass accretion rate 
was estimated to be ~4 orders of magnitude higher than 
the present value (Totani 2006). The current mass accre¬ 
tion rate of Sgr A* at the black hole horizon is Mbh ~ 
10~ 6 MEdd, here M^dd = lO^Edd/c 2 is the Eddington 
accretion rate (Yuan et al. 2003; see Yuan & Narayan 
2014 for a review of accretion models of Sgr A*). Thus 
the past accretion rate should be ~ 10 _ 2 MEdd- This 
rate is close to the largest accretion rate of a hot accre¬ 
tion flow, which is Merit ~ 0.07aMEdd ~ 2 X 10 _2 MEdd 
if a = 0.3 (Equation (27) in Yuan & Narayan (2014)). 
This indicates that Sgr A* was also in a hot accretion 
mode in the past. Specifically, for this accretion rate, 
the accretion flow consists of a thin disk outside of a 
transition radius r and a hot accretion flow inside of 
R tI . The value of i? t r is a function of accretion rate. For 
Mbh ~ 10 _2 MEdd) a reasonable value of R tr is ~ 200 R s 
(Yuan & Narayan 2004). This is the value we adopt in 
our fiducial model. 
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2.2. Winds from Hot Accretion Flows 

One of the most important findings in the theoretical 
study of hot accretion flow in the past few years is that 
strong winds exist (Yuan et al. 2012; Narayan et al. 2012; 
Li et al. 2013; Stone et al. 1999; Blandford & Begelman 
1999; Begelman 2012; Sadowski et al. 2013; see review 
by Yuan et al. 2015). This has been confirmed by ra¬ 
dio polarization (Aitken et al. 2000; Bower et al. 2003; 
Marrone et al. 2007) and X-ray observations (Wang et 
al. 2013). For example, it is found that only about 1% of 
the gas available at the Bondi radius finally falls onto the 
black hole while the rest is lost via winds. The detailed 
properties of winds, such as the velocity, mass flux, and 
angular distribution, have been obtained by magnetohy¬ 
drodynamic (MHD) simulations of hot accretion flows 
(Yuan et al. 2015; Yuan et al. 2012). Significant winds 
are found from 20P S up to the outer boundary of the 
hot accretion flow, which is P t r in the present case. The 
mass flux of winds is described by (Yuan et al. 2015): 


shown both analytically and numerically (Romanova et 
al. 1998; Blandford 2002; Hirose et al. 2004; Goodman 
& Uzdensky 2008; Uzdensky & Goodman 2008; Yuan et 
al. 2009). These events will efficiently accelerate parti¬ 
cles and produce CRs. Because winds are launched at 
the surface of the accretion flow (Yuan et al. 2015), we 
expect that many CRs must be contained in the winds. 
The ratio of the power of CRs (mainly CRp; the power of 
CRe can be neglected) and the total power of the winds 
(thermal gas and CRs) is defined as rjcn. — the “CR pa¬ 
rameter” . Since the details of such particle acceleration 
are still poorly understood, in the present work we treat 
rjcR as a free parameter. 

3. NUMERICAL SIMULATION 
3.1. Equations 

The MHD equations describing CR advection, diffu¬ 
sion, and dynamical coupling between the thermal gas 
and CRp are as follows, 


Af wind (P) = ^bh^t- (1) 

For Rtr ~ 200P S and Mbh ~ 0.02ME dd , we have 
M wind ~ 20%M Edd . We set M wind to be 18%Me<m in 
our simulations. The terminal poloidal velocity of winds 
is determined by the radius where the wind is launched 
(Yuan et al. 2015), 

«wind»(i2)«(0.2-0.4)t; fc (i2). (2) 

Here Vk{R) = (GM/R) 1 ' 2 is the Keplerian velocity at ra¬ 
dius R. At any given radius, say Rq, the poloidal velocity 
of winds is a mixture of winds launched from the region of 
r < Rq, and these winds have different poloidal velocities 
described by Equation (2). In the present work, we sim¬ 
ply assume that the poloidal velocity of winds is simply 
described by ~ Vk{Rt r) ~ 4.6%c since most of the mass 
flux of winds originate from around P t r- Thus the kinetic 
power of the thermal gas in the winds is 10 42 erg s -1 . The 
internal energy is negligible compared with the kinetic 
energy for the winds. 

2.3. Production of Cosmic Rays in the Accretion Flow 

The detailed acceleration process of cosmic rays (CRs) 
is beyond the scope of the present work. Here we as¬ 
sume that CRs are injected from the origin, together 
with the winds launched from the accretion flow of Sgr 
A*. These CRs may have the following origins. First, 
within the accretion flow the magnetic field is tangled, 
and the accretion flow is turbulent. Physical processes 
within the accretion flow such as MHD turbulence, mag¬ 
netic reconnection, and also weak shocks can accelerate 
some particles into relativistic energies. In fact, numer¬ 
ical simulations have shown the existence of magnetic 
reconnection in the hot accretion flows (e.g., Machida & 
Matsumoto 2003) and acceleration of particles has been 
studied (e.g., Ding et al. 2010). Second, at the coronal re¬ 
gion of the accretion flow, magnetic reconnection occurs 
perhaps more frequently. Loops of magnetic filed emerge 
into the corona from the accretion flow due to Parker 
instability. Since their foot points are anchored in the 
accretion flow which is differentially rotating and tur¬ 
bulence, reconnection occurs subsequently, as have been 


^ + pV • v = 0, (3) 

= -V(Pi + P 2 + Pb) - pV$ + V • T, (4) 
at 

—^ + V-(e 1 v) = -P 1 V-v + T:Vv + £ c , (5) 

at 

f)p 9 

— + V ■ (e 2 v - k ■ Ve 2 ) = -P 2 V • v, (6) 

HD 

— -Vx(vxB) = 0, (7) 

T = /n[Vv + (Vv) tr — ^IV • v]. (8) 

O 


Here Pi = (71 — l)ei and P2 = (72 — l)e2 are thermal 
pressure and CRp pressure respectively, in which 71 = 
5/3, 72 = 4/3 are adiabatic indices for ideal gas and 
CRp. Pb = B 2 /8n is the magnetic pressure and C c is 
cooling function. In order to avoid the density near the 
GC being too high, the gravitational force is adopted in 
the form of 


V<F = - 


2a 2 
r + r 0 


( 9 ) 


where r 0 is set to be 0.1 kpc. 

The diffusion coefficient tensor k can be written as: 


Kij 


s: _l Sij 


(«-L - K ||) 


BiBj 

B 2 ’ 


( 10 ) 


where k± and Ky are the diffusion coefficients perpendic¬ 
ular and parallel to the magnetic field, respectively. The 
relationship between the two coefficients is given by 


K± 1 

K|| l + (A||/r s ) 2 ’ 


( 11 ) 


in which Ay and r g are the parallel mean free path and 
the particle Larmor radius. In our case, Ay is much 
larger than r g , hence, <C km. Therefore we set k± 
to be 0 in most of our simulations. We should note that 
some works have shown that even if the turbulence level 
7 = (B 2 )/(P 2 + (B 2 )) (B is a random magnetic field) 
is small, say 0.1, the transverse diffusion coefficient k± 
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is much more effective than the quasi-linear result de¬ 
scribed by Equation (11), especially when the maximum 
scale of turbulence is comparable to the Larmor radius 
(e.g., Casse et al. 2001). In this case, the transverse diffu¬ 
sion can not be neglected. Considering the complexity of 
magnetic field coexisting with the CRs on the microscope 
scale, in one test we also include a non-zero transverse 
diffusion (refer to Table 1 for details). 

We include the cooling effect. We assume that the 
abundance of the halo is [Fe/H] = —0.5 (Miller & Breg- 
man 2013). The cooling function C c is set to be the 
cooling curve in collisional ionization equilibrium state 
(Sutherland & Dopita 1993). Cooling effect may be 
strong for central molecular zone (CMZ; Morris & Ser- 
abyn 1996, see paper I for its setup) gas, since the den¬ 
sity is very high there. But the heating effects for CMZ 
gas may also be strong, say radiative heating from stars 
and shock heating from supernovae, which are not con¬ 
sidered here. Here for simplicity, we assume that CMZ 
gas does not suffer from cooling and heating effects. We 
set a ground temperature of 10 4 K, below which cooling 
disappears. 

The viscosity is also included. As we have shown in 
paper I, if viscosity were not included in our model, both 
the Rayleigh-Taylor (RT) and Kelvin-Helmholtz (KH) 
instabilities would have developed into large-scale struc¬ 
tures during the formation of the Fermi bubbles. Viscos¬ 
ity may be affected by the magnetic field. Specifically, 
viscosity perpendicular to the magnetic field is strongly 
suppressed. However, we do not consider this complexity 
and still set an isotropic viscosity as in paper I and Guo 
et al. (2012). The dynamical viscosity coefficient /i is 
adopted to be 1.7 g cm -1 s _1 , which is slightly different 
from the value of 2 g cm -1 s -1 in the basic run in paper 
I. 

We neglect the influence of CRs on the magnetic field, 
which is complicated and difficult to be handled in our 
code. 

3.2. Simulation Setup 

We use the ZEUS3D code (Clarke 1996, 2010, also 
see http://www.ica.smu.ca/zeus3d/), to solve the above 
equations. ZEUS3D code is an Eulerian computational 
fluid dynamics code written in FORTRAN, and is a 
member of the ZEUS-code family. It can be used to 
deal with two-fluid problem. It adopts finite-difference 
method on a staggered grid with an accurate of the sec¬ 
ond order, and uses time-explicit, operator splitting so¬ 
lution in which the solution procedure is grouped into 
source step and transport step (see Stone & Norman 1992 
for details). For the diffusion process, we directly differ¬ 
ence the diffusion term in the source step. We adopt 
3-D Cartesian coordinates. Computational domain is al¬ 
most the same with paper I: —5.1 kpc ~ +5.1 kpc in 
X-, F-directions, and 0 kpc ~ +11.7 kpc in F-direction, 
and the X — Y plane is the Galactic plane while the Z- 
axis stretches along the Galactic pole. Sgr A* is just 
located in the origin. We adopt non-uniform grid, with 
Axi+i/ A Xi = 1.056 for x > 0, A%i/ A Xi = 0.9470 
for x < 0 , Ayj + i/ A yj = 1.056 for y > 0 Ayj + i/ A yj = 
0.9470 for y < 0, and Azk+i/ A Zk = 1.0285. The 
numbers of meshes are 120, 120 and 119 in X-, F-, in¬ 
direction, respectively. Since most of the parameters are 
the same with paper I, in the following subsections we 


mainly introduce the different parts. 

The global time step in the simulation is determined 
by the CFL condition (Clarke 1996). When viscosity, 
diffusion and cooling processes are considered, the time 
step is calculated by: 

dt = min{dt(CFL), dtvisc, dtdf, dtcooling} (12) 

in which the time step of viscosity process—dtvisc 
is calculated by p{ Ax) 2 / p (also see Guo et al. 
2012), diffusion process—dtdf is calculated by Cdf 
mm{dxidxjB 2 /(nuBiBj)} where Cdf is chosen to be 
0.25, and the cooling time step—dtcooling is neglected 
since it is much larger than the others. 

3.3. Initial Conditions 

The initial density distribution is set in the form of 
n e = rieoKpe (rk pc = r/1 kpc), where a is set to be 1.6. 
The values of n e o is set to be 7 x 10 -2 cm -3 except run C 
(see Table 1). The initial density is significantly higher 
than that in paper I, but close to the most recent work 
of Miller & Bregman (2015). Besides, we set the gradi¬ 
ent of the total pressure (thermal pressure plus magnetic 
pressure) to balance the gravitational force. 

Ptoti?) = Pi + Pb = 1.25 a 2 p(r) (13) 

where P tot , P\ and Pb (see section 3.4) are the total 
pressure, thermal pressure and magnetic pressure of the 
initial ISM. Therefore the initial ISM is not isothermal. 
We have also included the massive CMZ in our initial 
conditions (see paper I for the setup of CMZ), which 
plays an important role in collimating the winds. CRs 
are not included in the initial ISM, but injected from the 
GC (see Section 3.5). 

3.4. Magnetic Field 

We assume that the configuration of the magnetic field 
in GC in the past is the same with the present time. 
We also assume for simplicity that the winds launched 
from the accretion flow do not contain magnetic field. 
For the initial galactic magnetic field (GMF), we refer 
to the work of Jansson & Farrar (2012) for details. In 
our simulations, we assume the initial GMF contains two 
components, i.e., a large-scale regular field and a small- 
scale random field. 

The large-scale regular field can be divided into three 
components: a disk component, a toroidal halo compo¬ 
nent, and an out-of-plane component. The disk com¬ 
ponent is concentrated in the galactic disk, and drops 
into very weak magnetic field at height larger than 0.7 
kpc. Hence the disk component has a weak influence 
on the Fermi bubbles, while the toroidal and the out-of- 
plane component in the halo may play an important role. 
For simplicity, in our simulations, we only consider these 
two large-scale regular halo components, and superpose 
a tangled component as a small-scale random field: 

B = B tor + Bx + B t b ■ (14) 

where B tor is the toroidal halo component, Bx is the out- 
of-plane component, and B t b is the tangled field. The 
details of the three components are described in the Ap¬ 
pendix. 
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TABLE 1 

Simulations Parameters and Results 


Run 

n-eO 

(cm” 3 ) 

'Vwind 

(km/s) 

Mwind. b 

{M E dd) 

PCR c 

Kx/K|| 

^FB d 
(Myr) 

T e 

1 l 

(GeV cm -2 s 1 sr -1 ) 

L 7 * 

(erg/s) 

A 

7.3 X 10“ 2 

1.4 x 10 4 

18% 

50% 

0 

7.0 

3.3 x lO"® 

4.4 x 10 37 

B 

7.3 x 10 -2 

1.4 x 10 4 

18% 

50% 

0.01 

7.6 

3.8 x 10 -6 

5.0 x 10 37 

C 

3.7 x 10“ 2 

1.4 x 10 4 

18% 

50% 

0 

6.3 

1.3 x 10 -6 

1.7 x 10 37 

D 

7.3 x lO" 2 

1.4 x 10 4 

18% 

33% 

0 

7.9 

1.6 x 10 -6 

2.1 x 10 37 

E 

7.3 x 10“ 2 

7.0 x 10 3 

73% 

50% 

0 

7.6 

3.9 x 10 -6 

5.1 x 10 37 


a the velocity of winds. 

b the mass outflow rate of winds. 

c the ratio of the power of CR to the total power of winds (thermal gas + CR). 
d the age of the Fermi bubbles. 

e the gamma-ray intensity averaged on the projected map and integrated from O.lGeV to infinity. 
We have included the contribution from the secondary and primary leptons. 

f l he simulated gamma-ray luminosity of the Fermi bubbles at /f 7 >0.1 GeV. We have included 
the IC process as in / 7 . A fixed solid angle of 0.70 sr for both bubbles is assumed according to 
Ackermann et al. (2014). The center of each bubble is set to be ( R , \z\) = (0, 5 kpc). 


3.5. Cosmic Ray 

As we state in §2.3, we use a parameter rycR to describe 
the energy fraction of the CRp in the wind. We typically 
adopt TjcR = 50% in our simulations. But one lower 
value is also chosen in one model (see Table 1). 

The diffusion coefficient of CRs is energy dependent. 
A typical value by fitting CR data is Hi SO = (3 — 5) x 
10 28 cm 2 s -1 for ~1 GeV CR particles in a sufficiently 
tangled magnetic field on small scales, and scales as 
E 3 r 3 ~ E 3 r 6 , with E cr being the energy of a CR par¬ 
ticle (Strong et al. 2007). We simply treat the diffusion 
coefficients as independent of energy in our simulations. 
As argued in Yang et al. (2012), the diffusion coefficient 
Ki SO in a tangled magnetic field may be suppressed com¬ 
pared with K|| by a factor of 3. We therefore set k\\ to 
be 1 x 10 29 cm 2 s" 1 here, which is suitable for CRs at 
energies of a few GeV. 

4. RESULTS 

4.1. Distribution of CR protons 

The parameters of various models we have simulated 
are given in Table 1. We choose run A as our fiducial 
model, in which the age of 7 Myr is close to paper I. 
Figure 1 shows the distributions of density of thermal 
gas (n e ), energy density of CRp, and the temperature 
on X — Z plane. The density distribution here is similar 
to paper I. The density of the shocked ISM behind the 
forward shock is enhanced, and the density sharply de¬ 
creases crossing the contact discontinuity (CD; denoted 
by the purple dashed line in Figure 1) towards the inte¬ 
rior of the bubbles. Inside the CD, the bubble is filled 
with gas with very low density and very high temper¬ 
ature (a few 10 8 -10 9 K). The CRp pressure decreases 
away from the CD towards the exterior of the bubbles. 
Near the CD, the CRp pressure is comparable to the 
thermal pressure of the shocked ISM, and it expels the 
thermal gas away from the CD, leaving a zone with den¬ 
sity somewhat lower than the “typical” density of the 
shocked ISM. We call it a “permeated zone”. Therefore 
the shock structure in the present two-fluid MHD case is 
slightly different from that in paper I (see Figure 1 in pa¬ 
per I). The temperature in the shocked ISM ranges from 
3 x 10 6 K to 1 x 10 7 K when going from z = 3 kpc to 9 kpc. 
Since the radiation from the shocked ISM dominates the 


observed X-ray flux, we thus expect that observed X-ray 
spectrum should be fitted by the radiation from a multi¬ 
temperature medium. This seems to be consistent with 
the most recent Suzaku X-ray observation (Tahara et al. 
2015). 

The magnetic field in the shocked ISM is roughly par¬ 
allel to the CD, which agrees with the results of Yang et 
al. (2012). However, the alignment is not perfect, which 
is crucial for the hadronic scenario, as we will discuss 
later. Under such kind of magnetic field configuration, 
CRs can not diffuse too far away from the CD. Therefore 
the morphology of the bubbles is determined by the CD 
and this is also the reason why the edges of the Fermi 
bubbles look sharp. 

As shown by the middle plot of Figure 1, the distribu¬ 
tion of CRp looks like an “umbrella-like” structure. At 
the beginning, CRp move in the direction of the polar 
axis, and form the “handle” of the “umbrella”. After 
reaching the top, CRp fall down along the CD. This is 
because, CRp can only diffuse along the magnetic field 
lines which is roughly parallel to the CD. The density of 
CRp in the large volume enclosed by the umbrella-like 
structure is quite low, while the thermal pressure in this 
region is relatively high compared to the CR pressure 
(see Figure 2). This is of course because of the total 
pressure balance. A question is, the winds are blown 
out isotropically, why do CRp move like a “jet” at the 
beginning? Combining the middle panel in Figure 1 and 
Figure 2, we find that, the CRp “handle” along the polar 
axis is surrounded by a “wall” of high thermal pressure. 
The relatively high thermal pressure arises because of 
the strong viscous heating at the interface between the 
high-speed winds and the blown-up CMZ gas. In other 
words, the viscously heated high-pressure gas pushes the 
winds and form the “handle”. Besides, the magnetic field 
inside the winds (the purple region in the left panel of 
Figure 1) is parallel to the movement direction of the 
winds. This prevents the diffusion of CRs in the winds 
in the horizontal direction. 

We find that in the “umbrella face” the density of CRp 
decreases from the high to low latitudes, while it is on the 
opposite for the density of thermal gas. Our calculation 
shows that such a distribution can generate a roughly 
uniform gamma-ray emissivity shell in the sectional map 
thus a projection map with a constant brightness along 
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Fig. 1.— X — Z sectional views of the results of run A at t = 7 Myr. Coordinates are in units of kpc. Left : number density distribution 
of thermal electrons (n e ) in units of cm -3 . Velocity vectors are also plotted, with the color bar at the top of the plot denoting the value of 
velocity in units of km s — . Middle: energy density distribution of CRp (e2) in units of erg cm -3 . Magnetic field vectors are also plotted 
in this map in units of /iG. Right : temperature in units of Kelvin. The dashed lines in these three maps denotes the contact discontinuity 
(CD). 


the latitude (see §4.3 for details). If on the other hand 
the density of CRp did not increase with the latitude, 
the gamma-ray emissivity would have significantly de¬ 
creased with the height. This is because in either the 
hadronic or leptonic scenario, the density of both ther¬ 
mal gas and the radiation field decreases with the lat¬ 
itude, thus the surface brightness would have become 
dimmer with the increasing latitude. One advantage of 
this kind of “umbrella handle” structure is that it can 
keep the total energy of CRp from serious loss due to 
sideways adiabatic expansion. If on the other hand the 


CRp would follow a free and adiabatic expansion, more 
energy of CRp would have been lost, and it would be 
difficult to generate enough gamma-ray photons. 

As we have stated before, the magnetic field is not 
strictly parallel to the CD (see Figure 1). This allows 
some CRp to diffuse into the shocked ISM, form a per¬ 
meated zone outside the CD. We find that, about 60% 
of the total CRp have diffused outside the CD into the 
permeated zone, and the remaining 40% are still con¬ 
strained inside the CD. This is another crucial point for 
the hadronic scenario. The permeated zone is the domi- 
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Fig. 2.— Distribution of thermal energy density of run A at t = 7 
Myr. The values in the vertical color bar are in units of erg cm -3 . 
The horizontal color bar is for the velocity, which is denoted by 
arrows, and its value is in units of speed of light. 

nant region to generate the gamma-ray emission because 
of the relatively high density of thermal nuclei there, and 
only the CRp in the permeated zone are the “effective” 
energy source for the gamma-rays. 

4.2. Gamma-ray Radiation 

As we know, if the energy of the CRp is higher than a 
threshold value, the collisions between CRp and thermal 
protons can produce pions. The reaction channels are: 

p + p —> p + p + an 0 + b(n + + 7T _ ), (15) 

p + p —> p + n + n + + an 0 + b(n + + 7r“), (16) 

where a is generally equal to b and they denote the num¬ 
ber of pions produced in the reaction. Each neutral pion 
will instantly decay into two 7 -ray photons: 

7 r° -A 27 . (17) 

This process induces a lower limit of the gamma-ray pho¬ 
tons, which is about half of the energy of a neutral pion. 
Since the rest mass of a pion is 135 MeVc -2 , the gamma- 
ray spectrum will have a cutoff near ~ 70 MeV, which is a 
characteristic signature for the pion-decay, as observed in 
some supernova remnants (Ackermann et al. 2013). The 
gamma-ray spectrum of the Fermi bubbles also shows 
such a hardening feature below ~ 1 GeV, as highlighted 
in Crocker & Aharonian (2011). This is a possible evi¬ 
dence for the hadronic scenario. The charged pions pro¬ 
duced in pp collisions also contribute to the gamma-ray 
emission by generating high-energy secondary leptons. 
Charged pions decay in the follow ways: 

7 T + -» PL + + u p , /i + ->e + -|-[/ c + i( 1 , (18) 

7 t —> fi + p —> e + v e + G*- (19) 


All of these reactions can be finished instantaneously 
(the mean lifetime is 2.6 x 10^ 8 s for charged pions, and 
2.2 x 10 -6 s for muons). The secondary positrons and 
electrons are also of high energies, and they can scatter 
with the seed photons and produce gamma-rays. 

The products of pp collisions are calculated using the 
eparamlib package 5 , and the differential cross sections in 
pp collisions are from Kamae et al. (2006). We assume 
that the energy distribution of CRp follows such a power- 
law form: 

dn p (T p )/dT p = N PjPl T~ N (lGeV < T p < 2TeV), (20) 

where dn p is the number density of CRp with kinetic 
energies between T p and T p + dT p , N PtP i is a constant, T p 
is the kinetic energy of CRp, and the index N is set to be 
1.90. The number of final particles (including gamma- 
rays, secondary leptons and neutrinos) produced per unit 
volume, per unit time, and per unit energy in pp collision 
is given by (Ackermann et al. 2014): 


dQf 

lEf 


da(T p ,Ef) dn Pj T 

~ n HC dip , 
f dl P 


dE 


( 21 ) 


where the subscript / denotes the final particle species 
(such as 7 , e + , and e _ ), Ef is the energy of the final 
particle, da(T p , Ef)/dEf is the differential inclusive cross 
section (Kamae et al. 2006), hh is the number density of 
Hydrogen nuclei, and c is light speed. 

The secondary leptons undergo cooling via IC and syn¬ 
chrotron radiation. The evolution of the energy distribu¬ 
tion of secondary leptons can be obtained by: 


dn s {E s ,t) 

dt 


dE ., 


■[E s n s (E s ,t)] 


dQ s 

dE a ' 


( 22 ) 


where the subscript s means secondary positions or elec¬ 
trons, n s means the number density of the secondary 
leptons per unit energy interval, and E s is the energy 
of the secondary leptons. dQ s /dE s is the source term 
(see Equation (21)). The energy lost rate E s includes IC 
and synchrotron processes, with the magnetic field being 
2.5 pG according to Figure 1. We choose the stationary 
solution of this equation to calculate the IC scattering, 
and the seed photons are set to be the radiation field at 
(R,z) = (0, 5 kpc) (Porter & Strong 2005). 

As we have stated in §2.3, CRe can also be produced 
at the same processes as of producing CRp; thus there 
must be some CRe contained in the winds. The age of 
the Fermi bubbles in our model is comparable with the 
cooling timescale of CRe with energy of ~ 30-100 GeV 
(see Figure 28 in Su et al. 2010), so CRe with energy less 
than this value can still exist in the bubble. Based on 
these considerations, we can assume some primary CRe 
in the following power-law form: 

dn pe (E) / dE = N pe , pl E~ r e ~ E / E ™\ (23) 

where dn pe is the number density of the primary CRe 
with energies between E and E + dE : E is the energy 
of primary CRe, T is set to be 2.1, and the cut-off en¬ 
ergy E cut is set to be 60 GeV. The energy density of 
the primary CRe integrated from ~ 1 MeV to 60 GeV is 


5 https://github.com/niklask/cparamlib 
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7 x 10 -15 erg cm -3 , which is about 3 orders of magnitude 
lower than that of CRp. 

The number of gamma-ray photons produced per unit 
volume, time, and energy, which is the sum of tt° decays 
and IC scattering of the leptons, is given by: 


dNy 

dE-y 


dQ-y 

dEy 


+ IC (se + + se 


pe ), 


(24) 


IC(e) 



dcr jo (Ey , E e , Eph ) dn e driph 

—apjph 


dE~ 


dE e . 


dE, 


'ph 


(25) 


where se + , se~ and pe~ mean secondary positrons, 
secondary electrons and primary CR electrons, re¬ 
spectively. The cross sections of IC scattering, 
daic{Ey, E e , E p h)/ dEy is given by Blumenthal & Gould 
(1970), and we use the radiation field model at (i?, z ) = 
(0, 5 kpc) to set the seed photons dn p h/dE p h (Porter & 
Strong 2005). The total gamma-ray intensity, which is 
the gamma-ray energy observed by the telescope during 
per unit time, per unit energy in the vicinity of Ey, per 
unit solid angle, per unit area sensor is then 

I 7 = j : h dl = — jEy-^-dl, (26) 

where j 7 is the emissivity per unit solid angle, dl is 
the length element along the line-of-sight. As men¬ 
tioned above, / 7 contains three components: gamma-rays 
from 7r° decays I 7r o_ ) . 7 , IC scattering of secondary leptons 
J se _>. 7 , and IC scattering of primary CR electrons / pe _> 7 . 



0.1 1.0 10.0 100.0 1000.0 
E (GeV) 


Fig. 3. — The gamma-ray spectral energy distribution calculated 
based on run A. The rectangles with error bars show the latest 
observational results (Ackermann et al. 2014). The solid line is the 
sum of the dashed (for the n° decays), dotted (for IC process of the 
secondary leptons generated in hadronic reaction), and dot-dashed 
(for IC of the primary electrons) lines. 

From Equation (26), our calculation result of run A 
is shown in Figure 3, in which the vertical axis shows 
E 2 Iy. The result is averaged within the scope of the 
Fermi bubbles in the sky with latitudes larger than 10°. 
Note that we have multiplied the intensity of 7r° decays 
(dashed line) by a factor 1.5 as a rough correction for 
heavier ions (Mori 1997), i.e., the dashed line shows 
1.5-E 2 / 7r o_ s . 7 . Our model can fit the latest observation 


quite well. Specifically, the 7r° decays, the secondary lep¬ 
tons, and the primary CRe contribute 75%, 7% and 18% 
of the total intensity, respectively. 7r° decays dominant 
the origin of gamma-rays at Ey > 0.3 GeV; at Ey < 1 
GeV, IC of leptons is important. The leptons can signif¬ 
icantly soften the spectrum at Ey < 1 GeV, but still the 
hardening of spectrum in this band is clear, which can 
be regarded as the characteristic signature of 7r°-decay. 

The volume emissivity 47 tj 7 = EydNy/dEy for run A, 
B and E is shown in Figure 4. For simplicity, we only plot 
the emissivity of 7r° decays in this figure. The dashed line 
shows the CD. It is obvious that the gamma-ray emissiv¬ 
ity in the permeated zone outside the CD is much higher 
than that inside the CD. This is because the density of 
the shocked ISM is much higher outside the CD. The 
thickness of the permeated zone is typically only about 
0.6-0.7 kpc. Almost all of gamma-ray photons that are 
produced in n° decays come from this thin shell. There¬ 
fore, strictly speaking, it is the edge of the permeated 
zone caused by the CR diffusion, rather than the CD, 
that is the edge of the Fermi bubbles. 

4.3. The Gamma-Ray Surface Brightness 

The gamma-ray surface brightness is shown in Figure 
5. The contributions from both the 7r° decays and the IC 
(including primary and secondary leptons) are included, 
and the IC emissivity is assumed to be constant inside 
the bubble. The dashed and dot-dashed lines represent 
the edges of the observed north and south bubbles respec¬ 
tively (Su et al. 2010). The intensity is integrated from 
0.1 GeV to infinity, and is in units of GeV cm -2 s _1 sr -1 . 
The observation by Ackermann et al. (2014) masks the 
Galactic plane at |6| < 10°, so we also focus on |6| > 10°. 
The spatially averaged intensity for latitudes greater 
than 10° is 3.3 x 10 -6 GeV cm -2 s -1 sr -1 . From the 
figure we can find the following results. 

• The simulated morphology is close to the observa¬ 
tions, as in paper I. 

• The brightness does not vary with latitude, which 
is also consistent with observations. 

• The edge of the bubbles is not very sharp. Our sim¬ 
ulation shows that the angular width of the edge of 
the bubble is ~ 3°. This is consistent with the ob¬ 
servational result in Ackermann et al. (2014). This 
angular width is caused by the projection effect of 
the permeated zone with a typical thickness of 0.6- 
0.7 kpc. 

• We find that the surface brightness of the bub¬ 
bles is not very uniform, but shows a slight limb- 
brightened feature. The brightness near the edge is 
nearly twice the center (see the solid line in the bot¬ 
tom panel in Figure 5). This is because of the pro¬ 
jection effect of the permeated zone. In the work 
of Ackermann et al. (2014), the baseline model in 
1-3 GeV and 3-10 GeV bands also shows a blurry 
limb-brightened feature around b = 30° but with 
some uncertainties (refer to their Figures 22 and 
23). To compare with the observational data, here 
we calculate our gamma-ray flux in 3-10 GeV. The 
result is shown in the bottom panel of Figure 5 
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Fig. 4. — Simulated slices of the gamma-ray volume emissivity in E-y >0.1 GeV range. Here for simplicity, we only plot the contribution 
of 7T° decays. The three slices from left to right correspond to run A, B and E, respectively. The contour lines in each panels represent the 
density with the same meaning as in Figure 1. 



-30 -20 -10 0 10 20 30 

Longtitude 


Fig. 5.— Top : The gamma-ray intensity map for run A. Bottom : 
Flux averaged between b = 20° and 30° along the longitude. The 
solid and dashed lines are our simulation results, while the rect¬ 
angle with error bars and the shaded area are the observational 
results in 3-10 GeV band taken from Ackermann et al. (2014). 
The solid line is for E 7 > O.lGeV while the other lines and shaded 
area denote 3 times the flux. From the solid line, we can see the 
contrast between the bright edge and the dim center is about 1.6-2, 
and the boundary shows a width of about 3°. 


(dashed line). We can see that our result is roughly 
consistent with the data. Further observations are 
required to examine this issue. 


Yang et al. (2014) found that the morphology of the 
Fermi bubbles is energy-dependent: they are more ex¬ 
tended at high energies. From the figures in their paper, 
we estimate that the morphology at 10-30 GeV is about 


3°, or equivalently 0.5 kpc, larger than at 1 -2 GeV. Our 
model can explain this result. The physical reason is as 
follows. The above two gamma-ray bands roughly cor¬ 
respond to CRp with energy of ~ 200 GeV and ~ 20 
GeV, respectively. In our model, the total diffusion du¬ 
ration of CRp is about 7 million yrs (i.e., the age of the 
Fermi bubbles). Combing this duration, the respective 
diffusion speed of CRp at the above two energies, and 
the inclination angle of the magnetic field, we can calcu¬ 
late the diffusion depth of the CRs with the above two 
energies. We find that their difference is ~ 0.5kpc. 

4.4. The Effects of Changing Parameters 

In our fiducial model run A, we neglect the diffusion 
coefficient perpendicular to the magnetic field, i.e., n± = 
0. In run B, we include this transverse diffusion (see 
Table 1). In this case, we find that the permeated zone 
becomes slightly thicker (see the middle panel of Figure 
4), hence the width of the boundary of the bubble is 
slightly larger, and the contrast ratio of the bright edge 
and the dim center is smaller. The total luminosity is 
slightly enhanced (see Table 1). 

In run C, we have tested a case with the ISM density 
reduced by a factor of 2. We find that the gamma-ray 
intensity decreases by almost a factor of 3. 

Run D tests the effect of changing the “CR parameter” 
tjcr- Comparing runs A and D, we find that the final 
luminosity decreases by a factor of 2 when pcr decreases 
by one third. 

Increasing the mass flux in the winds but keeping all 
other parameters unchanged will obviously increase the 
gamma-ray intensity produced. Now we discuss the effect 
of another parameter, the transition radius between the 
hot and cold accretion flows Rt r . In our fiducial model, 
we set Rtr ~ 200I? S . If R tr is larger, the mass flux in the 
winds will be larger according to Equation (1), but the 
total kinetic power of thermal winds does not change ac¬ 
cording to Equation (2). In run E we increase the mass 
flux of the winds by a factor of 4 compared to run A while 
keeping the kinetic power of thermal gas in the winds un¬ 
changed. We find that, the gamma-ray luminosity only 
slightly increases. The increase is because more thermal 
nuclei are involved in the pp collisions inside the bubble. 
The increase is small since the gamma-rays are mainly 
produced in the permeated zone, whose properties are 
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not sensitive to the mass outflow rate but mainly to the 
power of winds. 

5. SUMMARY 

In a previous work (paper I) we have proposed an 
“accretion wind” model for the formation of the Fermi 
bubbles. Using hydrodynamic numerical simulations, we 
have shown in that work that the morphology of the 
bubbles can be successfully explained by the interac¬ 
tion between the winds launched from the hot accre¬ 
tion flow around the supermassive black hole Sgr A*. 
In the present work we continue this work by interpret¬ 
ing the gamma-ray emission of the bubbles. We invoke a 
hadronic model to explain the origin of gamma-ray radi¬ 
ation. In this scenario, the gamma-ray photons are pro¬ 
duced by the decay of neutral pions produced by the in¬ 
elastic collisions between CRp and thermal nuclei in the 
ISM. The CRp are produced within and perhaps more 
importantly at the surface of the hot accretion flow by 
processes like magnetic reconnection. These CRp will be 
carried by winds launched from the accretion flow and 
interact with the thermal particles in the ISM in the in¬ 
terface between the winds and ISM. The parameters of 
the winds are taken from the small scale MHD numeri¬ 
cal simulations of hot accretion flow (Yuan et al. 2015). 
The distribution of these CRp is calculated in the present 
work by three dimensional MHD simulations, considering 
the advection and diffusion of CRp in the thermal winds 
when the magnetic field is present. We then have calcu¬ 
lated the emitted gamma-ray spectrum from the bubbles 
based on the simulation data. We have compared the 
results with observations and found that our model can 
not only explain the morphology of the bubbles but also 
explain the gamma-ray spectrum very well (Figure 3). 

Some details of the results can be summarized below. 

• Although the winds are almost isotropic, the CRp 
do not expand isotropically into the galactic halo. 
Rather, an “umbrella-like” structure is formed 
(Figure 2). The CRp first move along the “handle” 
of the “umbrella”, after they reach the top of the 
bubbles, they diffuse downwards along the surface 
of the “umbrella”. In this scenario, CRp occupy 
a relatively small volume, thus the total energy of 
CRp does not lose much. This kind of distribu¬ 
tion is also conducive to form a constant surface 
brightness in the latitude direction. 

• We find that the magnetic field is slightly mis¬ 
aligned with the CD, which makes most of the 
CRp diffuse across the CD into some depth of the 
shocked ISM, forming a permeated zone (Figure 1). 
The thickness of the permeated zone is only about 
0.6-0.7 kpc, but is the dominant region for gener¬ 
ating the gamma-ray photons since thermal nuclei 
are much denser than inside the CD. Therefore it 
is the edge of the permeated zone, rather than the 
CD, that defines the edge of the Fermi bubbles. 

• Our model can fit the observed gamma-ray spec¬ 
trum quite well (Figure 3). The gamma-ray mainly 
comes from the decay of neutral pions produced in 
the pp collisions. 

• Our accretion wind model can explain the width 
of the boundary of the Fermi bubbles, which is 


caused by the projection effect of the permeated 
zone (Figure 5). 

• Our result shows a somewhat limb-brightened sur¬ 
face brightness, with the contrast between the 
bright edge and the dim center of 1.6-2. This is 
roughly consistent with the observed blurry limb- 
brightened feature in 1-3 GeV and 3-10 GeV maps 
in Ackermann et al. (2014) (bottom panel of Figure 
5). 

6. DISCUSSION 

6.1. The Giant Radio Loop Structure 

In paper I we have mentioned that our model can ex¬ 
plain the ROSAT X-ray features. Now we discuss the 
interpretation of another structure, the giant radio loop 
called Loop I. This structure was first observed in the 
radio band (Large et al. 1962), and most recently also 
detected in gamma-ray (Ackermann et al. 2014). Loop I 
radio image roughly surrounds the north Fermi bubble, 
and the recent gamma-ray observation also shows that 
the gamma-ray image of Loop I happens to surround 
both the north Fermi bubble and the lower part of the 
south Fermi bubble (refer to Figure 13 in Ackermann et 
al. 2014). The gamma-ray spectral index is a ss —2.4, 
which is significantly softer than that of the Fermi bub¬ 
bles (Su et al. 2010, Ackermann et al. 2014). The origin 
of Loop I is still not clear. 

The North Polar Spur (NPS) which was observed in 
both the radio and X-ray bands, is the brightest region 
in Loop I. Guo & Mathews (2012) propose that this 
structure shares the same origin with the Fermi bub¬ 
bles and corresponds to the shocked ISM. Here we focus 
on the gamma-ray emission of Loop I (or NPS). The 
forward shock outside the Fermi bubbles can acceler¬ 
ate both electrons and protons. The CRe could emit 
radio emission through synchrotron radiation. Gamma- 
ray emission could be produced through IC scattering 
of soft seed photons by CRe, or pion decays through pp 
collisions. Quantitatively, if the gamma ray emission is 
dominated by IC of CRe, the observed gamma-ray spec¬ 
tral index requires that the energy spectral index of CRe 
is roughly p ss —2.4. If these electrons are mainly ac¬ 
celerated through Fermi first-order acceleration in the 
forward shock, the required Mach number of the shock is 
~ 4 (Drury 1983). The Mach number in our simulation 
is ~ 4 — 7, close to the required value. 

6.2. Winds versus Jets 

Since the accretion flow around Sgr A* has been a hot 
accretion flow for millions of years, a jet must also exist 
together with the winds (Yuan & Narayan 2014). In our 
“accretion wind” model of the Fermi bubbles, we do not 
include jet. We assume that the role of jet is negligible 
to the formation of the bubbles. If we do not consider 
complexity such as the precession or turbulence of ISM, 
the jet will simply pierce through the ISM in a narrow 
low-density channel through which the jet can freely flow, 
without producing a bubble-like structure (Vernaleo & 
Reynolds 2006). 

On the other hand, in the literature a type of jet model 
has been proposed to explain the formation of the Fermi 
bubbles (e.g., Guo & Mathews 2012; Guo et al. 2012; 
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Yang et al. 2012, 2013). Several assumptions are adopted 
in this model. Firstly, to explain the bubbles which is 
perpendicular to the Galactic plane, the jet is required 
to be perpendicular to the Galactic plane as well. This 
assumption seems to be too strong since radio surveys 
show this is generally not the case (Zubovas et al. 2011). 
Secondly, the required mass lost rate in the jet is too 
large. It is ~ 0.3MEdd in the hydrodynamical model 
(Guo et al. 2012) and ~ 300-MEdd in the MHD model 
(Yang et al. 2012, 2013), respectively. Theoretical models 
of jet formation predict that only a small fraction of the 
accretion matter can go into the jet, say ~ 10% for the 
“disk jet” or much smaller for the “Blandford-Znajeck” 
jet (see Yuan & Narayan 2014 for the review of various jet 
models and Figure 5 in Yuan et al. 2015 for the mass flux 
in the disk jet). As we state in §2.1, the accretion rate in 
the underlying accretion flow in the past was likely only 
~ 10 _2 MEdd- Thirdly, the velocity of the jet required in 
the Yang et al. (2012, 2013) is very low, only ~ 0.03c. In 
Guo et al. (2012) it is ~ 0.1c, which is still much smaller 
than the numerical simulation and some observational 
results, unless significant deceleration of jet occurs. For 
example, the speed of the jet in M87 is > 0.986c (e.g., 
Biretta et al. 1999). 

Our result has interesting implications to the formation 
of cavities and bubbles often observed in galaxy clusters. 
The general explanation is that they are formed because 
of the interaction between jets and ISM. The usual argu¬ 
ment is that we have observed a jet there. However, we 
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APPENDIX 


In this part we introduce the three components of the GMF adopted in our simulation. 

Toroidal halo component: Bt or (see Jansson & Farrar 2012). This component only has one direction: </> in the north 
halo, or -<j> in the south halo. The strength of the field is in the form: 

Btor = Btor 0^ I ^ ° ^ Hdisk disk) iX B(v, T n , IC/i)), (1) 

L(z,h,w) = ( l + e - 2 (l*l-'0/«')- 1 , (2) 

in which zo = 5.3 kpc, hdisk = 0.4 kpc, Wdisk = 0.27 kpc, r n = 9.2 kpc and Wh = 0.2 kpc. We set the B tor o = 1.4 pG 
which is the value of northern halo. Considering that we have not included the disk magnetic field component, here we 
set L(z, hdisk, Wdisk) to be a constant of 1, which is a very good approximation for z > 0.7 kpc, otherwise the toroidal 
magnetic field would be too weak near the galactic plane. 

Out-of-plane component: B x (see Jansson & Farrar 2012). The field lines of the out-of-plane component are straight 
lines in the space, and show a mirror symmetry about the Z = 0 plane (see the middle panel in Figure 6 for this 
magnetic field). This magnetic field lines look like “X-shaped”, and is similar to the “X-shaped” field lines in some 
edge-on galaxies, such as NGC 891, NGC 4666, NGC 5775 and so on (Haverkorn & Heesen 2012). If we set the 
cylindrical radius where the out-of-plane field line intersecting the plane of Z = 0 to be R = r pi the magnetic field 
strength on the Z = 0 plane is 

bx{r p ) = B X0 e~ r »' rx , (3) 

where Bxo is 4.4 /rG, and r x is 2.9 kpc. The space can be divided into two parts, divided by a critical cylindrical radius 
R = r c x on Z = 0 plane and a constant elevation angle of 49° starting from this critical circle. The elevation angles 
are different in the two parts. Specially, for r p > r x , the elevation angle is a constant of 49° , which is represented by 
9 X . The magnitude of out-of-plane magnetic field in r p > r\ region is described by: 

\Bx\ = b x (r p ) r -f , 


p tan 6 X ’ 

where r c x is a constant of 4.8 kpc. 

For r p < r x , the elevation angle linearly changes from at R = r v to 0 at R = 0. The magnetic field is described 
in the form: 

\B x \=bx(r p )C-^) 2 , 
r 

r ■ r x 

rp ~ 

I z\ 

6\(r , z) = arctan(—-—) 
r-r p 

Tangled field: B t b■ The direction of the magnetic field B t b is set to be a spatially periodic form to mimic the 
“tangled” status, in the form of: 


( 6 ) 

( 7 ) 

( 8 ) 


( 4 ) 

( 5 ) 


B\ = B 0 [sm(K 2 y + A 2 ) + sin(A' 3 ,z + A 3 )], 
B 2 = -Bo[sin(-ffix + Ai) + sin (K 3 z + A 3 )], 
B 3 = B 0 [sin(iFix + A 3 ) + sin [K 2 y + A 2 )], 


( 9 ) 

( 10 ) 

(11) 
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where Bq is a constant of 0.5 /rG, -?G, 2,3 are set to be 6.28 kpc -1 which means that the initial GMF is “tangled” on 
length scale of 1.0 kpc in each direction, and the phases Ai 2,3 are set to be 0.5, 0.5 and 1.0, respectively. To explore 
the influence of the phases, we also try another set of data: Ai ^.3 = 1.5, 1.5, and 1.0, respectively. In this case, the 
magnetic held on the galactic polar axis is relatively stronger, but we find that the final result is almost the same. It’s 
easy to prove that magnetic held with this form satishes the zero divergence. 

The three components are shown in Figure 6. The large-scale regular magnetic held is the sum of B tor and Bx ■ 



Fig. 6. — Magnetic field configurations on different slices. Left: purely large-scale regular magnetic field (without the “tangled” component) 
on the slice of Z = 4 kpc. Middle: purely large-scale regular magnetic field (without the “tangled” component) on the slice of Y = 0 kpc. 
Right: purely tangled magnetic field on the slice of Z = 4 kpc. Color bars show the strength of magnetic field in units of //G. 

































































